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We discuss results of numerical analyses of the 3+1 dimensional Nambu-Jona-Lasinio 
model at non-zero baryon density. By studying relevant order parameters, it is shown that 
the model goes through a transition between a vacuum with broken chiral symmetry, and 
a chirally restored phase at high baryon chemical potential p with a diquark condensate 
roughly proportional to the area of the Fermi surface and a broken U(l) baryon number 
symmetry. The condensation of diquark pairs from the fermionic spectrum is shown to lead 
to an energy gap A about the Fermi energy Ef, which is found to be w 15% of the vacuum 
fermion mass and roughly independent of \i in the chirally restored phase. The ground-state 
is believed, therefore, to be that of a BCS superfluid. We also discuss the effect of simulating 
on a finite volume, and any caveats this places on the conclusions made above. 

§1. Introduction 

The physics of compact stars is a field that has attracted much interest in recent 
years. At zero temperature, deconfined quark matter within the core of these objects 
is believed to be unstable with respect to a ground-state in which quark (and hole) 
pairs in the region of the Fermi surface condense out of the spectrum in a manner 
analogous to that in the BCS theory of superconductivity; because diquark pairs 
cannot be colour singlet, this phenomenon is known as colour superconductivity 
(CSC)P The simplest CSC ground-state, known as 2SC, is one in which condensa- 
tion occurs between the approximately degenerate u and d quarks, leaving an energy 
gap 2A about the Fermi surface in their spectra. The fact that the force between 
quarks in QCD is strongly attractive means that estimates for the magnitude of A 
can be as large as 50-100MeV, i.e. large enough to remain present at typical compact 
star temperatures of 0(lMeV)P^ 

The confirmation of this picture by a first principles lattice QCD calculation, 
however, remains elusive due to the difficulties of performing Monte-Carlo simula- 
tions with baryon chemical potential fi ^ 0. Whilst the "sign problem" remains, 
the only way of probing the cold, dense region of the QCD phase diagram non- 
perturbatively is by using model field theories, such as 2 colour QCD or the Nambu- 
Jona-Lasinio (NJL) model. 

In both of these models, the functional weight remains both real and positive, 
even with fi ^ 0®® suc h that they are simulable using standard lattice techniques. 

*' This work was completed with the support of PPARC and in collaboration with Simon Hands. 
Thanks also go to Tom Kingaby for enlightening discussions. 
**' Current address. 
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Lattice simulations of the NJL model show that chiral symmetry restoration occurs 
at some critical chemical potential \i c ~ Uq, where Uq is the constituent quark mass, 
in qualitative agreement with analytic approaches such as the large- iV c (or Hartree) 
approximation.® Furthermore, whilst this model has no gauge degrees of freedom 
and cannot describe confined matter, the fundamental quarks that play the role of 
the baryons do obey the correct statistics, such that the effect of increasing fi beyond 
He is to build up a Fermi surface. This should be contrasted with 2 colour QCD, 
in which the qq baryons are bosonic and the onset of matter occurs at /j, c ~ m n /2. 
Together, these features make the lattice NJL model an ideal tool with which to study 
Fermi surface phenomena, such as CSC, at phenomenologically relevant densities. 

§2. The model 

The model studied herein is the 3+1 dimensional lattice NJL model, which in 
continuum notation has the Euclidean Lagrangian 
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where m is the bare- or current-quark mass, and the source j has been added to 
allow the measurement of a baryon number violating diquark condensate in a finite 
volume system. 

When formulated on a discrete lattice of sites x with separation a between 
nearest neighbours, the model is described in terms of the staggered fermion fields 
X, X, C an d C by the action 
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where A is the Nambu-Gor'kov matrix 
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In turn, M is the fermion kinetic matrix 
1 
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xy 2a 



(x,x) 



(2-4) 



where T pq is a vector of the Pauli matrices which act on the Nf = 2 isospin degrees 
of freedom. The lattice model also exhibits an additional N c = 8 "colour" degrees of 
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freedom, due to staggered fermion doubling and the choice of (j2-2[) . The auxiliary 
scalar and pseudoscalar fields a and 7? are defined on sites x of a dual lattice, (x, x) 
represent the sum over the 16 dual sites neighbouring x, and the symbols Tj u (x) and 
e(x) are the phases (— l) xoH H^-i and (— i^o+^i+^+zs respectively. 

In the limit that m and j — > 0, (|2-2j) has the same global symmetries as QCD 
with two massless flavours, i.e. the SU(2)l <8> SU(2)r chiral symmetry of left- and 
right-handed quarks 

X -» (PlU + P R V) X ; X^x(PlV^ + PrU^); 
( ^ (P L V + P R U)C; (^((PlUI + PrVI); (2-5) 
(er + Z7r.f) — > V(<r + ijf.r)U^, 

where Pl,r = \ (1 ± e ( :r )) and C7 and V are elements of SU(2), and the ?7(1)b 
symmetry of baryon number 

X -» e iQ x, X -» e ~ ia X- (2-6) 

As in QCD, at fi = chiral condensation leads the spontaneous generation of a 
constituent quark mass Eq 3> m which breaks SU(2)l®SU(2)r down to the SU(2)j 
of isospin. These reassuring features suggest that by studying the model's global 
symmetries and their breaking, we can learn about the analogous symmetry breaking 
patterns in QCD. 

In 3+1 dimensions the NJL model is only an effective theory, such that the lattice 
parameters must be chosen to match low energy phenomenology. By applying the 
large- N c (Hartree) approximation analytically to staggered quarks, we determine a 
suitable set, leading to physically reasonable results, to be 
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§3. Order parameters 



In previous work,® we have performed simulations of the model outlined in SJ2 
in an attempt to identify diquark condensation in the dense phase. In particular, we 
measure the diquark condensate 

1 dlnZ 1 /„ f t 2 



where Z is the partition function, and compare it to the chiral condensate 

M'v-toT-wV'K-u V)- (3 ' 2> 



We perform simulations for a range of \x on V = Lt x lattices with Lt = L s = 12, 
16 and 20 using a standard hybrid Monte-Carlo (HMC) algorithm, in which the 
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source strength is set to zero during the update of the auxiliary background fields, 
but allowed to vary over a range of j = j* 7^ during the quenched measurement 
of the quantities above. An interesting empirical observation is that results in the 
chiral sector are found to scale linearly with inverse volume, whilst diquark observ- 
ables scale with inverse temporal extent, which in Euclidean space corresponds to 
the temperature. Accordingly, (qq) is extrapolated to V^ 1 — > 0, whilst (qq) is ex- 
trapolated first to L^ 1 — > and then j € [0.3, 1.0] — ► 0.*) The results are plotted in 
Fig. ^ As expected, at \i = chiral symmetry is broken leading to a non-zero chiral 

08 I 1 1 1 1 1 

0.7 - (qq) I — 9 — I 

0.6 - Area of Fermi surface - 




Fig. 1. Order parameters as functions of fj,. The solid line represents the area of the Fermi surface 
calculated in the large- N c limit. 

condensate whilst {qq) ~ 0. As chemical potential is increased, the system goes 
through a crossover into a phase with approximate chiral symmetry restoration and 
a non-zero diquark condensate. In particular, (qq) is found to increase roughly as the 
area of the lattice Fermi surface, plotted in arbitrary units for free fermions in the 
large- N c limit as the solid curve. This supports the naive picture of the condensate 
being a measure of the density of pair states about the Fermi surface, contributing 
to diquark condensation. 

§4. The superfluid gap 

In the previous section, we presented evidence for superfluidity in the form a 
non-zero diquark condensate in the high-/x phase. This evidence must be treated as 
indirect, however, since this quantity is not measurable in an experiment. Further- 
more, even if the sign problem were overcome, one could not measure (qq) in lattice 
QCD, since in this theory the colour superconducting phase is characterised by the 
breaking of a gauge symmetry; Elitzur's theorem states that one cannot define a 
local order parameter to distinguish the existence of such a phase in a gauge invari- 

*' Although all diquark observables are measured for 0.1 < j < 1.0, the behaviour of data with 
j < 0.3 appears to disagree markedly with that of those at higher j, and these are ignored for the 
purposes of these fits. We attribute this to finite volume errors, which we discuss briefly in 
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ant way.^ Instead, we turn to more direct evidence in the form of the global order 
parameter for the BCS phase, the energy gap A. 

We do this by studying the fermionic dispersion relation E(k), extracted from 
the Gor'kov time-slice propagator 

G(k, t) = J2 -4 _1 (0, 0; x, t)e~ iLs , (4-1) 

on L t x L x x Ly z lattices with L t = 16, 20 and 24, L x = 96 and L VyZ = 12. The 
momenta sampled are k = (2Tm/L x ,0,0) for n = 0, 1, 2, . . ., implying that for our 
choice of lattices there are 25 independent momentum modes between k = and 
7r/2 in the k x direction. We measure Q at jx = 0.8 for a range of j, using standard 
lattice methods for meson correlators, in the same partially quenched approximation 
employed in the measurement of the order parameter. It is found empirically that 
whilst in general Q has 16 complex components in Nambu-Gor'kov space, only two 
parts of these are independent and non-zero. These are 

N(k,t) = Re(G 13 (k,t)) k A(k,t)=Im(g 21 (k,t)), (4-2) 

which in the limit that j — > correspond to the propagators for particles 

N(k,t) ~ (q(x)q(y)) n (4-3) 

and for superpositions of particles and holes 

A(k,t)~( q (x)q(y)) 21 . (4-4) 

We label these the "normal" and "anomalous" propagators respectively. Fitting 
these to 

N(k, t) = Ae~ Et + Be- E(Lt - t} (4-5a) 

A(k, t) = C(e~ Et - e E{Lt - ty ) (4-5b) 



we extract the amplitudes A(k,j), 
B(k,j) and C(k,j) and the excitation 
energy E(k,j), and as with (qq) ex- 
trapolate to L^T 1 — > and then j € 
[0.3,1.0] -> 0. 

Fig. [21 shows the amplitudes A, B 
and C as functions of j in the limit 
T, j — > 0. At low k, A^> B correspond- 
ing to a predominantly forward moving 
signal, which implies that the normal 
propagator deep within the Fermi sea is 
dominated by hole-like excitations. At 
large k the reverse is true, implying that 
excitations above the Fermi surface are 




k 

Fig. 2. Propagator amplitudes at fi = 0.8. 
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dominated by particles. The point at which A{k) and B(k) coincide, therefore, can 
be identified with the momentum at the Fermi surface, /cpP» Of particular interest is 
the fact that the coefficient for the anomalous propagator C(k) is non-zero in a broad 
peak about kp, even in the j ' — » limit, which implies that there are particle-hole 
mixed states with indefinite baryon number, an indirect signal for superfluidity via 
a BCS mechanism. 




k k 
Fig. 3. Lattice dispersion relation and typical free fermion dispersion relation at /i = 0.8. 



The left-hand panel of Fig. |31 shows the fermion dispersion relation at fi = 0.8, i.e. 
the energy of the fermionic excitations as a function of their momenta, compared to 
that of non-interacting lattice fermions. The free fermion dispersion relation clearly 
has two branches, representing hole excitations for k < kp and particle excitations 
for k > kp. In between one observes that £'(sin _1 (sinh(/ia))) = 0. In contrast 
to this, in the dispersion relation extracted from the lattice data one can see no 
discontinuity between the two branches, another sign of particle-hole mixing. More 
importantly, at no point does E{k) pass through E = and there is a distinct gap 
between this point and the minimum. 

This can be seen in a more familiar light if one plots the hole branch as nega- 
tive, which is done in the right-hand panel of Fig. [3 This makes the free fermion 
dispersion relation a smooth continuous curve, which is similar to those observed in 
lattice theories with no BCS gapP® whilst for our smooth dispersion relation this 
introduces a discontinuity at k ~ kp signalling an energy gap 2 A. This is the first 
direct observation of a BCS gap in a lattice simulation. 

In order to learn something of the chemical potential dependence of the gap 
we perform a naive extrapolation of data extracted from Lt = 16 and 20 lattices 
to zero temperature and assign a conservative estimate of the error; the Lt = 24 
data required for a full statistical treatment such as that applied at fj, = 0.8 are too 
expensive to reproduce for a range of [i in the chirally restored phase. Data are 
then extrapolated to j — > as before so that we can plot the dispersion relations 
and gain an estimate for the gap. These are plotted in the left-hand panel of Fig. |1J 
It is worth noting that whilst a linear extrapolation through two points is of little 
statistical value, the results produced at \x = 0.8 with and without using the Lt = 24 



A BCS Gapped Superfluid on the Lattice 



7 




data are found to be consistent for all k. 

The estimates for A(fj,) are plotted in the right-hand panel of Fig. which 
shows that A appears roughly constant for fj, > fi c , whilst {qq} rises with the area 
of the Fermi surface. This is consistent with the simple-minded assumption that the 
order parameter counts the density of states from within a shell of thickness O(A) 
about the Fermi surface that contribute to diquark condensation, such that in the 
continuum 

(qq) oc Ail 2 . (4-6) 

Finally, in order to express our best measurement of the gap, Aa(fia = 0.8) = 
0.053(6), in a form independent of the lattice spacing a, we extract the vacuum 
fermion mass Uq from simulations with Lt = 16, 20 and 24 and find 

A = 0.15(2). (4-7) 

Assuming that Uq ~ 400MeV, this implies that in physical units, the size of the gap 
is A ~ 60MeV, which is consistent with the predictions of Ref. |2J) • 



§5. Finite volume effects 



An important point that we have failed, thus far, to highlight is that the con- 
clusions of the previous sections rely on the disregarding of data with j < 0.3, for 
which the order parameter and other diquark observables deviate sharply from the 
fitted form. We attribute this to finite volume errors, which are particularly large 
in this system j*» a fact we believe to be due to the difficulty of representing a small 
shell of states close to the Fermi surface on a discrete momentum-mode lattice. To 
address this issue, parallel simulations are currently being attempted on very large 
lattices, with the aim of showing that the on such systems the deficiency at low-j 
is reduced. Thusfar, these simulations have proven hard, with very small time-steps 



8 



David N. Walters 



and trajectories required for acceptance rates of over 50%. This might be interpreted 
as a sign of encouragement, however, as the failure of the algorithm may be a signal 
that the partially quenched approximation is insufficient for such large volumes. If 
the volume is large enough that the proposed finite size effects are overcome, there 
may be enough momentum modes within the shell about the Fermi surface for di- 
quark condensation to occur even with j = 0, as partial quenching enforces during 
the field updates. The algorithm would then be attempting to represent a system 
with an exact Goldstone mode and encounter the standard finite size effects related 
to divergent Goldstone fluctuations. A full simulation, with j ^ during all steps 
of the algorithm, may yet be required to resolve this issue. 

§6. Summary 

To summarise, our evidence for s-wave superfluidity via a BCS instability is of 
some importance, since this is the first time that the presence of such a phase has 
been demonstrated in a relativistic quantum field theory using a systematic calcula- 
tional technique. Although the 3 + Id NJL model is only a simplistic effective field 
theory, this work can be interpreted phenomenologically as non-perturbative evi- 
dence for BCS colour superconductivity in QCD with two degenerate flavours, since 
these two theories have the same global symmetry structure. This is of particular 
importance whilst the persistence of the sign problem prevents the numerical solu- 
tion of full SU(3) QCD. Also, the fact that our results agree with the predictions of 
self-consistent treatments of this simple model adds credence to solutions of similar 
models with more complicated flavour structures and interactions. 

In future work it would be interesting to further this analysis to study the stabil- 
ity of the superfluid phase by performing simulations either at non-zero temperature 
with the aim of measuring the critical temperature of the superfluid phase, or with 
the Fermi surfaces for "up" and "down" quarks separated via the introduction of 
a small non-zero isospin chemical potential. The latter option may be difficult, 
however, since initial studies indicate that with this introduction, detM becomes 
complex. 

The most pressing issue, however, is to resolve the finite volume effects discussed 
above, since it is only once this issue has been addressed that the conclusions of J|l 
& ?2]may be trusted completely. 
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